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ABSTRACT 


Two important aspects of Hg/Cd/Te crystal growth processes are discussed herein. First, 
the thermal field and, second the fluid movement in the melt zone. The thermal analysis 
includes numerical calculation of axisymmetric heat conduction within the sample. It 
also includes a three-dimensional radiation model to calculate the radiative heat 
exchange between the furnace and the crystal as determined by the complex geometry of 
the furnace and the adiabatic shield. The thermal analysis also includes a crystal 
conductivity which is dependent on temperature and composition. To tackle the fluid 
flow aspect of the problem, an attempt was made to use a newly developed 
incompressible flow code based on the slight compressibility, and hence the finite sound 
speed, of all real fluids. 
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INTRODUCTION 

This contract (NAS8-36483) started on April 24, 1985 with an intended period of 
performance of 12 months. During these first 12 months, a successful thermal analysis 
of Hg/Cd/Te crystal growth was performed. The analysis consisted of numerical 
axisymmetric heat conduction within the sample. A three-dimensional radiation model 
was incorporated to model the radiative heat exchange between the furnace and the 
crystal as determined by the complex geometry of the furnace and the adiabatic shield. 
The interface between the solid and melt phases was considered to be a non-isothermal 
surface determined by data from the phase diagram of the Hg/Cd/Te mixture. The 
conductivity in both the melt and solid phases was calculated as a function of 
temperature and composition. The objective of this phase was to determine the effect of 
furnace temperatures on the shape and location of the interface, and hence the quality 
and composition of the grown solid crystal. The effect of a metal tube enclosure on the 
thermal field within the crystal was also investigated. The results of this first phase of 
the study were completely documented and delivered to NASA as an interim report on 
May 1, 1986. A black and white copy of this report is included in Appendix "A” for the 
convenience of the reader. 

As a result of the satisfactory completion of the first phase, the contract was extended 
to December 1986 to initiate a new phase of the study. This new phase is concerned with 
the fluid movement in the melt zone and its effect on the quality of the solid crystals. 

Continuum, Inc. developed an incompressible flow code for use in a separate contract 
(NAS8-35508) as well as the second phase of this effort. The code incorporates a slight 
compressibility approach which takes advantage of the bulk modulus and finite sound 
speed of all real fluids. The finite element numerical analog uses a dynamic differencing 
scheme based, in part, on a variational principle for computational fluid dynamics. 
Details of the code's theory and calculation procedure are given in the following section. 

The development and testing of the new code occupied most of the extended period of 
performance. An attempt was then made to use the code for calculating melt zone 
circulation. The attempt was unsuccessful because of two main reasons, namely, the 
extremely small magnitude of velocity compared to pressure, and the excessive amount 
of computational time required by the code. In other words, the approach is suitable for 
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relatively high speed flow regimes and requires expensive main-frame computer 
resources which were not available under this contract. The code was delivered to NASA 
under Contract NAS8-35508. Details of the theory and calculation procedure are given 
in the following sections, as well as conclusions and recommendations for future work. 


- 2 - 



icr 1 0* 


CI-FR-0102 


INCOMPRESSIBLE FLOW CODE DEVELOPMENT 


Nomenclature 


a 

c 

E 

E B 

F 

I 

P 


s 

t 

T 

u,v,w 

U 

V 

6 

Vi 

u 

5 

P 

T 


-speed of sound 
-incompressible specific heat 
-total system energy 
bulk modulus of elasticity 
-flux integrals, Eq. (14) 

-identity tensor 
-pressure 
-velocity vector 
-heat transfer 
-entropy 
-time 

-temperature 

x-, y-, z- components of velocity, respectively 
-conserved variable array, Eq. (14) 

-volume 

artificial compressibility 
viscosity 

kinematic viscosity 
allocation parameters, Eq. (16) 
density 
stress tensor 
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Literature Review 


The solution of the governing equations in primitive variable form for incompressible 
flows has proven to be a difficult task. The number of different equation sets used is 
nearly as large as the number of researchers working in the field. This is in sharp 
contrast to the solution of compressible gas flow in which, with the ideal gas assumption, 
the equations are well known. For incompressible gas and/or liquid flows, there is no 
comparable equation of state with which to close the system of equations. Density, and 
consequently pressure, are usually eliminated from the continuity equation, which 
becomes a constraint on the momentum equations. Although the pressure does appear in 
the momentum equations as spatial derivatives, neither of these equations can be solved 
for pressure without introducing directional bias (Ref. 1); nor can a transient type of 
solution be developed. Thus, from the outset one is faced with "N" equations for "N + 1” 
unknowns. 

Early numerical solutions to incompressible flow avoided this problem by using the 
vorticity-stream function formulation, which is thoroughly treated in Ref. 1. This 
procedure suffers from many well known limitations, particularly for internal and three- 
dimensional problems (see, e.g., Refs. 1,2). With the realistic three-dimensional 
problems which are within reach by using current and future generations of 
supercomputers, the use of these equations will be severely limited. 

One of the first procedures for obtaining pressures of incompressible flows required the 
solution of a Poisson equation (Ref. 1). This equation has the form 

v 2 P = S (1) 


where S is a non-linear function of the velocities and their first three derivatives. Also, 
Neumann boundary conditions are required to solve Eq. (1). The solution to Eq. (1) must 
be obtained iteratively at each step, a very time consuming process. Another drawback 
is that the solution does not behave well when complex geometries are being studied. 
Furthermore, the solution of Eq. (1) only once at each time step does not properly 
account for the velocity-pressure coupling, as demonstrated in Ref. 3. Since it is well 
known, Eq. (1) was specifically studied for inclusion into the crystal growth model; 
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however, because of the problems mentioned above, it was determined to be a non-viable 
alternative. 

An increasingly popular solution used in finite difference applications is the artificial 
compressibility developed by Chorin (Ref. 4). This formulation assumes an equation of 
state as 

P = p/6 (2) 


where 6 is the artificial compressibility. Using Eq. (2) the continuity equation becomes 
s( — V v • q = o (3) 

Eq. (3) allows for an explicit finite difference solution for P; however, since S is 
technically a relaxation parameter and has no physical basis, unsteady solutions are 
meaningless. Steady state solutions, on the other hand, have been obtained successfully 
by numerous researchers. Chang and Kwak (Ref. 5) and Kwak, et al (Ref. 6) present one 
the of the most in depth and most successful applications of the method, including 
applications to SSME problems (Ref. 7). 

A Finite Element Method (FEM) which embodies the same general principle as artificial 
compressibility is the penalty function formuation (Ref. 8). This procedure is widely used 
in pure FEM techniques, although many variants have been introduced in an attempt to 
obtain accurate results (for example, Refs. 9, 10). This further serves to illustrate the 
difficulty of obtaining pressures for incompressible flows. 

The artificial compressibility method was considered for use in developing the crystal 
growth model, however, the essentially arbitrary nature of 6 in Eq. (3) creates several 
problems. Perhaps the major difficulty is in the use of the variational scheme explained 
later in this report. With no physical basis for S , the differencing procedure did not 
behave welL Other problems include the trial and error procedure required to determine 
6 , and the unavailability of an unsteady solution. 
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There are numerous other methods for treating the velocity-pressure coupling, such as 
the transformed equations of Solomon and Szymczak (Ref. 11) or the assumed elemental 
"deviatoric stress-velocity-pressure" method of Yang and Atluri (Ref. 12). Few of these 
methods provide a set of physically accurate equations for the primitive variables. Also, 
transient problems can, at best, only be "studied” by iterating for successive steady state 
solutions, which is very costly and inefficient. Numerical damping is also necessary to 
obtain solutions in most of the available methods, including the artificial compressibility 
method. 

A recently developed procedure which attempts to overcome most of these deficiencies 
is presented in Ref. 13. This work develops the pressure solution based on a finite sound 
speed, in other words, a compressible formulation. The compressibility is used to obtain 
an equation of state, which in turn introduces the pressure into the compressible form of 
the continuity equation. The resulting equation is 

3 p 2 - - 

— + q • V P + pc (V • q) = 0 (4) 

at 

Eq. (4) has been successfully applied in Ref. 13 to both transient and high Reynolds 
number problems. 

The procedure to be used in developing the crystal growth model is very similar to that 
of Ref. 13. The final equation differs somewhat from Eq. (4) because fewer assumptions 
will be made. The details of the analysis are presented below. 
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Basic Equations 

In the development of the compressible formulation, the integral form of the equations 
of motion are used. Several advantages of this formulation occur, as pointed out in Ref. 
14. The fundamental equations are: 


continuity: 



(5) 


momentum: 


cv 9 1 


y/*[(pq)q + 

CS 


(Ip- T )]*ds = 0 


( 6 ) 


energy: 


cv g t 



(IP-t) 


q + Q ] • ds = 0 


(7) 


For this analysis, as with all incompressible methods, the assumption is made that P 4 P 
(T), hence the energy equation (Eq. 7) decouples from the other equations. Since liquid 
flows are to be studied, it should also be pointed out that gravity, if assumed to be 
important, is easily included in Eq. ( 6 ) as a source term. 


- 7 - 



CI-FR-0102 


Pressure Solution 


As previously discussed, the current approach to the pressure solution is to take 
advantage of the finite sound speed in all fluids, liquid or gas. The procedure is more 
correctly called a slight compressibility method, rather than an incompressible method. 

An equation of state is assumed to exist in the form 

P = P(p) (8) 


The compressiblity of any fluid is expressed by the bulk modulus of elasticity, 


dP 



P 


(9) 


and the finite acoustic speed in a liquid is given by 


2 

a = 


e b dP 


The equation of state may now be written as 
dP = a 2 d p 


( 10 ) 


(ID 


To derive the pressure equation, Eq. (11) is primarily used in differential form, however, 
a reference state must also be assumed for integration purposes. This reference state, 
p r and p r , must be assumed to be the usual theoretical incompressible values of P and 
p for the fluid being considered. Substitution of Eq. (11) into Eq. (5) now yields the 
pressure equation. 
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fff— dV + ^Pq^i + ( p a 2 - P ) ffq-4s = 0 (12) 

CV at cs cs 

The convective part of Eq. (12) is written as two terms to facilitate discussion of some 
salient features of the method. First, it should be noted that as the sound speed, a , 
approaches infinity (the true incompressible case) Eq. (12) reduces to the standard 
incompressible form of continuity. Another feature is that the artificial compressibility 
method can be viewed as a limiting case of Eq. (12). If the contribution of the pressure 
to the convective terms is discarded and p r a 2 = 1/6 , then Eq. (12) reduces to the 
integral form of Eq. (3). Perhaps the most important feature is the physical reasoning 
used to obtain Eq. (12). In addition to more correctly describing the physics of the 
flowfield, this method allows the use of a variationally based dynamic differencing 
scheme. 

Eq. (11) could (and perhaps should) be included in the momentum equations. The error 
introduced by using p = p r in Eq. (6) is essentially zero for steady state solutions, 
which can be viewed as incompressible solutions. For unsteady solutions, the error is of 
order ^P/pa 2 vs. p and, since a is large, this error is assumed to be negligible. 


Numerical Analog 

The numerical integration used to solve Eqs. (6) and (12) is based in part on the 
variational procedure developed by Prozan in Refs. 14-16. This procedure will be 
summarized below, and the application to the current work will then be discussed. 


The compressible transport equations may be generalized as 


^/// 

3 1 cv 


U 


n 


dV = - F 

n 


j n — 1 , • • . , 5 


(13) 


where 


U = 


P-) 
pq J 

►Tj 

II 

( pq • ds _ 

J [pqq + (ip - t) ] 

pE 

CS 

( 'pEq + (IP - t ) 


q + Q ] • di 


(14) 
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Now, the domain of integration is subdivided into E finite elements and Eqs. (13) become 
E E 

I UU n 4V/ A t) e = - I F (15) 

e= 1 e=l ’ 


A functional distribution of the variables is assumed over the surfaces of the elements 
and the flux terms F n 0 are evaluated using values from the previous time step. 


The (AU n AV/At ) g represent the total accumulation of the nth conserved quantity 

within element e. To obtain the conserved variables at a new time step, nodes are placed 

at the corners of each element and the (AU AV/At) are allocated to the various 

n e 

nodes. Thus, the conserved quantities at the nodes are determined by assembling the 
contributions from its surrounding elements. This accumulation can be written 


(AU n AV/At). 


k 

r r . F 

g =1 ^nje,] n , e 


(16) 


where j is the node number, k is the total number of elements surrounding node j, and the 

£ . are the allocation parameters. 

s n,e,j 

The assembled equation, Eq. (16), is analogous to a finite difference expression in which 

the spatial transformations are numerically embedded in the analog. It may be 

interpreted as a general form of the finite difference scheme since different finite 

difference algorithms can be derived with selected allocation parameters. The 

variational differencing procedure differs from other schemes in that it does not dictate 

a fixed r . throughout the course of integration, but rather, changes the 
n , e , j 

parameters dynamically in both time and space according to the variational principle. 
This principle requires that the rate of entropy production be maximized. The 
variational differencing scheme thus uses the transport equations as equality constraints 
and dynamically determines the allocation parameters to achieve maximum stability of 
the system. 

To apply the numerical analog to the current investigation, a procedure for determining 
the allocation parameters similar to the procedure of the Ref. 14 must be developed. A 
first attempt of this development was presented in Ref. 17. This method was not 
entirely satisfactory, and further research was performed. The entropy functional of 
Ref. 17 was further manipulated and the following derivatives were obtained: 
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3pS 1 1 q 2 

= —rr clnT - c + - 

9P a 2 T 


where c is the specific heat, 
dps u 

3pu T 

3ps v 

3pv T 

3 ps w 

dpw T 


(17) 


(18a) 


(18b) 


(18c) 


Using Eqs. (17) and (18), the allocation parameters of Eq. (16) are determined in the same 
manner presented in Ref. 14. This procedure provided good stability for many cases. 
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RESULTS 


The results of the first phase of this study were completely documented and delivered to 
NASA as an interim report on May 1, 1986. A black and white copy of this report is 
included in Appendix "A" for the convenience of the reader. In the second phase of the 
project an attempt was made to use the newly developed incompressible code for 
calculating melt zone circulation. The attempt was unsuccessful for two main reasons. 
First, the velocity magnitude is so small that the numerical noise dominates the 
flowfield. To illustrate, consider the typical term (P + pu 2 ) which models the effect 
of pressure and inertial forces in the momentum equation. The velocity u is in the order 
of microns per second, while the pressure is in the order of 10® dyne/cm 2 . Thus, the 
round-off error in the pressure term may be larger than the magnitude of the cross- 
velocity term. It is therefore, impossible to calculate the velocity field without 
resorting to special procedures. 

The second difficulty is the excessive amount of computational time required by the 
code. Since a realistic speed of sound in liquids (5000 ft/sec) must be used, the time step 
size required for stability is very small. The funds available for this project were not 
sufficient to secure any amount of main-frame computer resources to perform the 
required calculation. 

The following conclusions have been arrived at: 

1. The molten column is stable since it is heated from above and since the heavier 
material exists at the bottom of the column next to the interface. 

2. The melt movement is due only to gravity which causes the liquid to "slide” along 
the curved interface. 

3. The temperature and pressures inside the sample are much higher than the 
velocity magnitude. 

4. Calculation of the velocity field in the molten crystal requires the development of 
a specialized code using an implicit rather than an explicit approach. 
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RECOMMENDATIONS 


It has been concluded that a strong effort should be devoted to developing a procedure 
for calculating the molten movement and its effect on the quality of the crystal. The 
following is a list of recommendations and guidelines for such an effort: 


1. A "small perterbation" form of the governing equations must be used. 

2. A thorough review of literature should be made to identify computer codes 
suitable for this kind of problems. 

3. Effects of gravity, magnetic fields and density variations must be accurately 
accounted for. 

4. An implicit algorithm must be developed to solve the governing equations. 
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ABSTRACT 


A thermal analysis of Hg/Cd/Te crystal growth in a Bridgman Cell is made using 
Continuum’s VAST code and workstation facilities. The analysis consists of numerical 
axisymmetric heat -conduction within the sample. A three-dimensional radiation model is 
incorporated to model the radiative heat exchange between the furnace and the crystal 
as determined by the complex geometry of the furnace and the adiabatic shield. The 
interface between the solid and melt phases is a non-isothermal surface determined by 
data from the phase diagram of the Hg/Cd/Te mixture. The conductivity in both the 
melt and solid phases is a function of temperature and composition. The objective of this 
study is to determine the effect of furnace temperatures on the shape and location of the 
interface, and hence the quality and composition of the grown solid crystal. The effect 
of a metal tube enclosure on the thermal field within the crystal is also investigated. 
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INTRODUCTION 


Mercury Cadmium Telluride crystals are essential for manufacturing infrared detection 
devices. The purity and composition uniformity of the crystal are most important and, 
unfortunately, difficult to achieve. The objective of this work is to identify methods and 
optimize conditions leading to the production of uniform and pure crystals. 

This investigation specifically addresses the thermal phenomena which dominates the 
process of crystal growth. Numerical axisymmetric heat conduction within the 
Hg/Cd/Te sample is performed. The sample is regarded as consisting of a solid zone and 
a molten zone. The thermal conductivity in both zones is temperature and composition 
dependent. The interface between the two zones is taken in the first phase as the 706 °C 
isotherm, and in the second phase of the research is defined as a non-isothermal surface 
whose temperature distribution is obtained from the Hg/Cd/Te phase diagram. The 
composition in the solid zone is considered uniform if the interface is flat. If the 
interface is curved, the solid crystal composition is non-uniform in the radial direction 
and is taken from empirical data supplied by NASA. The heat exchange between the 
furnace and the sample, through the air gap and the glass ampule, is a complex 
phenomena of radiation and conduction. This exchanged heat is accurately calculated 
herein through a three-dimensional determination of the view factors, and a heat balance 
on the glass wall elements. Thus, the thermal energy flux at the sample’s surface is 
determined and used as a boundary condition for solving the heat conduction equation 
within the sample’s domain. The solution provides the thermal field within the sample 
and hence the shape of the interface and the composition distribution in the solid 
crystal. Parametric cases were calculated to study the effect of the furnace operating 
conditions and geometry on the quality of the crystals. 

The following two sections of this report summarize the thermal properties and data used 
throughout the study, and the details of the boundary conditions calculation. Next, 
results of parametric cases (assuming isothermal interface) are presented. Finally, more 
advanced calculations of the interface's shape and the solid composition are presented 
followed by results of more parametric cases. At the end of the report, a list of 
conclusions and recommendations, a list of references and an appendix of mathematical 
details are attached. 
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THERMAL PROPERTIES 


The Hg-^_ x Cd x Te sample is assumed to have a constant density of 7.45 gm/cm 3 . Density 
variations and convective motions in the molten zone are not considered in this study. 
The initial bulk composition, x Q is 0.2 for all the parametric cases presented herein. The 
density of the glass amplule is 2.203 gm/cm 3 and its specific heat is 1046.7 x 10 6 
mm^/sec 3 °K. The thermal conductivity of the ampule’s glass, kg» its radiative 
transmissivity, x, and the air's thermal conductivity, k ft , are all functions of temperature 
as shown in Fig. 1. The Stefan-Boltzmann constant is 5.6696 x 10 -3 gm/sec 3 °K^, and the 
emissivities of all furnace surfaces are equal to 1. 


The local composition, x, is calculated at each node in the liquid phase using the 
following one-dimensional species solution (1): 


x = x„ 
o 


S-l 

1 - exp(-RZ/D) 

S 


(for Z > 0) 


where S is a segregation coefficient equal to 3.9; R is the growth rate (assumed equal to 
the translation rate); D is an effective diffusion coefficient equal to 5.5 x 10 -5 cm 3 /sec; 
and Z is the vertical distance above the interface whose location is defined later in this 
report. According to the above equation, x is equal to x Q /S at Z = 0 and changes 
exponentially in Z to approach x Q as Z becomes large. In the solid phase (Z < 0), the 
composition is first assumed constant and equal to x 0 « In the last stage of this 
investigation, the solid composition is radially varied in a manner dependent on the shape 
of the interface. Details are given elsewhere in this report. 


The thermal conductivity, K, is calculated at each node using the following equations 
(supplied by NASA): 


K = a p C p 

a = A q - A x T* + A 2 T* 2 - A 3 T* 3 (solid phase) 
a = B q ln(T) - Bj (liquid phase) 
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Fig. 1 Properties of Air and Glass at Different Temperatures 
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where ot is the thermal diffussivity in mm^/sec; Cp is the heat capacity of the Hg Cd Te 
mixture (constant at 1.8428 x 10^ cm^/sec^ °K); T is the temperature in °C; and T* is 
T/1000. The empirical constants, A 0 , Aj, A 2 , A 3 , B Q , and Bj are given in Table 1 as 
functions of composition. 


The pseudobinary phase diagram of Hg Cd Te mixtures, Fig. 2, is used herein for 
interface and solid composition determination as explained elsewhere in this report. The 
diagram is obtained from Ref. (2) and is described by the following empirical equations: 


s 


X 


s 


X 


Cj s i n ( 


" * 
- T ) 
2 


71 

+ C« s i n ( — 
1 2 


+ c 4 ][7 

x® ( T - 670 ) / 20 



) + logiQ (9 T 
(for T > 690 °C) 
(for T < 690 °C) 


+ 1) 


= D, T 


+ D 2 T 


*2 


+ D 3 T 


*3 


D. T 
4 


*4 


where x s is the solidus composition; x 1 is the liquidus composition; T is temperature in °C; 
and T* is equal to (T - 670)/412. The constants Cj and Dj (i = 1,....4) are given in Table 2. 
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X 

A 0 

Al 

A 2 

a 3 

Bo 

Bl 

0.000 

2.784 

9.938 

20.98 

16.21 

10.4655 

67.401 

0.052 

1.719 

6.065 

14.48 

8.797 

7.25493 

46.767 

0.107 

1.512 

6.494 

14.81 

11.03 

7.09724 

46.024 

0.205 

1.125 

4.568 

11.03 

8.427 

6.44051 

42.041 

0.301 

1.031 

4.652 

11.59 

8.933 

5.69691 

37.423 


Table 1 Values of the Empirical Conductivity Coefficients (NASA) 


i 

C i 

Di 

1 

0.502804 

0.607640 

2 

0.165390 

0.077209 

3 

0.746318 

0.696167 

4 

-0.413546 

-0.381683 


Table 2 Values of the Phase Diagram's Empirical Constants (Ref. 2) 
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THERMAL BOUNDARY CONDITIONS 


Thermal energy is exchanged between the alloy and the furnace, through the glass wall of 
the ampule, by radiation and conduction. A part of the radiated energy is directly trans- 
mitted to the alloy while the glass wall absorbs the remaining part and reflects a negli- 
gible amount. The radiant energy absorbed by the glass, along with the energy conducted 
through the air gap, is conducted through the two-dimensional domain of the glass wall. 
A part of this conducted energy reaches the alloy, while the remaining part flows in the 
longitudinal direction. The net result is a certain amount of energy reaching or leaving 
the outer surface of the alloy per unit area per unit time. An accurate evaluation of this 
energy flux is vital for predicting a correct temperature distribution within the alloy. 
Once calculated, this energy flux, Q, can easily be used to impose the proper boundary 
condition on the solution of the energy equation. 

In order to calculate the boundary heat flux, Q, a "lumped” heat balance is performed on 
the portion of the glass wall associated with each boundary node (see Fig. 3). This 
element of glass is assumed to have a uniform temperature T fl which is an average of the 
node temperature, T n , and the temperature of the outer surface of the ampule at a point 
facing the node, T T n is obtained from the solution of the energy equation at the 
previous time iteration and Tg is calculated by equating the heat conducted to and 
through the air/glass interface: 


(K /5 ) (T,- T ) 
a a f g 


(K /5 ) (T - T ) 
g g g n 7 


where and Kg are the conductivities of air and glass; <5 a and 6 g are widths of the air 
gap and glass wall; and Tf is the furnace temperature at a point facing the node. In a 
similar manner, T Qa and T ab are calculated to represent average temperatures in the 
glass portions above and below the portion under consideration. 
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The set of temperatures Tf, Tg, T a , T n , T aa and T ab is used to calculate the 
following heat fluxes: 


Q = (K /6 ) (T - T ) 

^a g za' aa a' 


«b - ( V { zb> (T ab- V 


Q =(2K/6 ) (T-T) 
v g g g Si 


V (2 W (T n" V 


Q r (absorbed) = Q r ( 1 - t) 


Q (transmitted) = Qt 
r r 


where t is the transmissivity of the glass, and Q r is the radiative heat flux exchanged 
between the furnace and the glass element. Details for calculating Q r are given in 
Appendix A. In each of the above formulae, K g is the glass conductivity at a 
temperature equal to the average of the pair of temperatures used to calculate the 
flux. Also, t is the transmissivity of the glass at a temperature equal to (T^ + T a )/2. 


The above energy fluxes cause the temperature of the glass element to change from T a 
to T fl ». The new T fl ’ is calculated from the heat balance on the glass "lump": 


(V g P g C g /At) (T a * T a ) = (Q a + V A 1 + (Q g + Q r (absorbed) )Ag+ 

where Vg is the volume of the glass element; P is the density of the glass; Cg is the 
specific heat of glass; and A t is the time step size. The areas Aj_, A2 and A3 are shown 
in Fig. 3. Finally, the flux Q which represents the energy exchange between the furnace 
and the outer surface of the alloy in the vicinity of the node is calculated as: 

Q = (2Kg/6g) (T^ - T n ) + Q r ( t r ansmi t ted ) 


CI-IR-0087 


At each boundary node, the energy per unit volume P E 0 ]<j (obtained by solving the 
energy equation) is adjusted to accommodate the flux Q as follows: 

pE „ew * pE old + Q M At/v 

where v is the volume of the alloy element associated with the node. This energy 
adjustment is performed at the end of each time step. 
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INTERMEDIATE RESULTS 


The heat conduction equation is solved in the axisymmetric domain descretized as shown 
in Fig. 4. The thermal boundary condition at the wall nodes is as described in the 
preceding section. The temperature at the top end nodes is fixed at the hot furnace 
value, and the temperature of the bottom end nodes is equal to the cold furnace 
temperature. The interface is defined as the 706 °C isotherm. Thus, the molten and the 
solid zones are. identified and separated by a sharp interface. The solid phase 
composition is constant at 0.2. Molten phase composition, and the conductivity of both 
phases are determined as described earlier. 

Under the above conditions, 15 parametric cases are calculated. The first three cases 
are intended to show the effect of the difference between the hot and cold furnace 
temperatures (THOT-TCOLD) on the shape and distribution of the steady-state isotherms 
within the sample. It is easily seen in Fig. 5 that the sample's isotherms become closer 
to each other as the difference between THOT and TCOLD increases. In other words, a 
steeper temperature gradient across the heat barrier results in a steeper temperature 
gradient in the sample around the interface. The translation rate in all the three cases is 
0.15 lJ/sec. 

Runs 4 through 11 are designed to investigate the effect of TCOLD on the location and 
hence, the shape of the interface. It is seen in Fig. 6 that the interface moves upward as 
TCOLD decreases while THOT is held constant. When the interface is in the lowest 
position (Run 4) it curves downward to become concave relative to the melt zone. As it 
moves upward, the interface flattens while TCOLD is in the range (425, 300 °C) and then 
curves upward to become convex relative to the melt zone as seen in Run 11. A closer 
look at the interfaces of Runs 6, 7 and 9 reveals, however, a subtle change in the 
interface curvature for the TCOLD range of (425, 300 °C). As seen in Fig. 7, the 
interface - while shifting steadily upwards as TCOLD falls - changes from concave to 
flat and then to concave again before it eventually becomes convex as in Runs 10 and 
11. It is believed that the interplay between radiation, conduction and conductivity 
changes is responsible for this behavior of the interface. A detailed conceptual 
explanation is found in Ref. 1. 
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Fig. 5 Effect of CTHOT - TCOLD) on Thermal Gradients in the Sample 
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Fig. 7 A Closer Look at Interface Shape 
(Runs 6, 7 and 9) 
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Runs 12 through 15 show the effect of the sample’s translation rate on the steady-state 
isotherms and the shape and location of the interface. As seen in Fig. 8, the translation 
rate must be changed dramatically in order to cause a notable change in the distribution 
of the sample’s isotherms. 
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ADVANCED CALCULATIONS 


O 

In this stage of the research, the interface is no longer considered as the 706 C 
isotherm. Also, the solid phase composition is no longer fixed at the x Q value. The 
calculations of the melt phase composition and the thermal conductivity of both phases 
remain as explained earlier in this report. 

Reference 3 presents a comprehensive empirical study of the solid crystal composition. 
Based on this study, it is herein assumed that the composition in the solid phase varies 
radially (and not axially) within the range (0.19, 0.235) for an initial bulk composition of 
0.2. From the two composition values, 0.19 and 0.235, the temperature range within 
which solidification must have occured is easily obtained from the phase diagram as 
shown in Fig. 9 (711.5, 702.5 °C). This temperature range is divided into ten segments 
which is equal to the number of elements used to descretize the sample in the radial 
direction (see Fig. 4). The set of temperatures T^, T 2 , ..., T^j is then used to find the 
location of the interface on the axial planes passing through the corresponding nodes. 
For example, on the axial grid line number 6, the temperature of the interface is Tg. 
Hence, the location of the interface is easily determined by finding Tg on this axial grid 
line. Repeating the process over all the 11 grid lines, the interface is completely 
determined at the end of each time step. 

If the interface is concave, as in Fig. 10(a), a disk of melt moving downwards will solidify 
first at the axial grid line number 1 at a temperature Tj. Thus, for the next time step, 
the interface is located by finding Tj on line number 1; T 2 on line 2 and so on until Tjj 
on line 11 which is the center line. In this case, the solid composition is assumed to vary 
from Xj (Fig. 9) at grid line number 1 to xjj at grid line 11. 

If the interface is convex, as in Fig. 10(b), the molten disk moving downwards will 
solidify first at the axial grid line number 11 under a temperature of Tj. Thus, for the 
next time step, the interface is located by finding Tj on line number 11; T 2 on line 
number 10 and so on until T^ on line number 1. In this case, the solid composition 
varies from Xj at grid line number 11 to xjj at grid line 1. 
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The interface is flat when the set of temperatures Tj, Tjj all lie on a horizontal 

line within certain tolerance. In this case the solid composition is uniform at 0.2. 
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FINAL RESULTS 


Six parametric cases (Runs 16 through 21) are calculated under the conditions mentioned 
in the preceding section. In runs 16, 17 and 18, THOT is 812 °C and TCOLD is 610, 350, 
and 75 °C respectively. The translation rate in all the three runs is 0.15 p/sec. Figure 
11 shows the effect of TCOLD on the shape of the sample's isotherms and on the 
composition distribution in both the molten and solid phases. Figure 12 is a detailed view 
of the isotherms, in the solidification zone with the non-isotherm al interface shown as a 
dashed line. It is seen that if THOT is fixed at 812 C, a TCOLD of 350 °C is required to 
produce a flat interface and, hence, a uniform composition in the solid crystal. A 
numerical optimization tool for Hg Cd Te crystal quality control is thus created. Runs 
similar to 16, 17 and 18 can be easily used to predict the quality of the grown crystals 
under a variety of furnace operating conditions and geometries. 

Runs 19, 20 and 21 differ from the three preceding runs in that the glass ampule is 
enclosed in a metal tube of 12.94 mm diameter. The disk-like heat barrier surrounds this 
metal tube on the outside and ensures that the temperature distribution along the tube is 
identical to that of the furnace. Thus, the procedure for calculating the radiative heat 
flux at the sample's boundary remains the same. Only the furnace is now at the position 
of the metal tube and the heat barrier is absent. 

Figure 13 shows the results of Runs 19, 20, 21 compared to those of Runs 16, 17, 18. It is 
seen that the effect of the metal tube is to bring the isotherms at the sample's boundary 
closer together. That is expected since, in effect, the furnace was moved closer to the 
sample forcing its temperature distribution to follow more closely that of the furnace. 
Note also that in order to obtain a flat interface, a TCOLD of 400 °C is needed instead of 
350 °C. 
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Fig. 11 (Cont.) Sample’s Isotherms and Composition 

Under Different Cold Furnace Temperatures 
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Fig. 13 Effect of Metal Enclosure on Sample’s Isotherms 
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CONCLUSIONS AND RECOMMENDATIONS 


An accurate thermal analyzer of Hg Cd Te crystal growth was developed which includes 
the following aspects: 

a. Axisym metric heat conduction in the sample. 

b. Three-dimensional calculation to model the radiative heat exchange between the 
furnace and the crystal as determined by the complex geometry of the furnace 
and the adiabatic shield. 

c. Non-isothermal interface between the solid and molten phases determined by the 
phase diagram of the Hg Cd Te mixture. 

d. Temperature and composition dependent conductivity in both the melt and solid 
regions. 

The model was used to investigate the effect of the furnace temperature distribution, 
and the translation rate on the location and shape of the interface and, hence, on the 
quality of the grown crystal as indicated by its composition distribution. Also, the effect 
of a metal tube enclosure was investigated. The results are displayed in Figs. 5 through 
13 and discussed earlier in this report. The model can be systematically used to predict 
the quality of the crystal grown under different conditions. 

It is recommended that a new phase of the study be started to add the following 
attributes to the model: 

a. The fluid mechanics of the motion of the dense liquid in the molten zone; its 
effect on the shape of the interface and the quality of the grown crystal. 

b. Effect of gravity (or its absence) on the fluid motion. 

c. Effect of any magnetic field that may exist and affect the flow in the melt zone. 
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APPENDIX A 


BASIC CALCULATION FOR RADIATIVE EXCHANGE 



The basic expression for the radiative exchange between the furnace and a certain 
location within the glass or the alloy’s surface can be written as follows: 


q dA f -dA n = £ f 




-dA_ 


dA f - 


e a T 
n n 


d V dA f 


dA 


(A-l) 


This is the net energy diffusely leaving the differential area dAf and reaching the 
differential area dA n per unit time (see figure A-l). It is assumed that the media 
between dAf and dA n does not participate in the radiation process. A negative quantity 
means that the net energy leaves dA n and reaches dAf . o is the Stefan-Boltzmann 
constant; (£f, Tf ) and (e n , T n ) are emissivities and temperatures of dAf and 
dA n respectively. The view factors between the two differential areas are given by: 


F dA f -dA n = C0S *f C0S<t> n dA n 
F dA n -dA f = C0S *f C0S< ^n dA f 


/ TT S 2 , 
/ IT S 2 , 


(A-2) 


( A-3 ) 


where S is the distance between the two differential areas; 4>f and 4> n represent the 
angles between the line S and the normals to dAf and dA n respectively. If 
(R n , © n , z n ) and (Rf, 0f, z f ^ are the cylindrical coordinates of the centers of 
dA n and dAf respectively, the following expressions are easily obtainable: 


dAf = R f d0 f dz f , (A-4) 

dA_ = R d0 dz , ( A-5 ) 

n n n n 7 

S 2 = R 2 - 2 R f R n cos ( © f - e n ) + R n + (z f‘ z n )2 » (A-6) 


A-l 




cos<j) f = (1/S) [R f - R n cos(0 f - e n )] , 


( A-7 ) 


cos<J> n = (1/S) [R f cos(6 f - 6 n ) - R n ] . 


( A-8 ) 


Using Equations (A-2) through (A-8), Equation (A-l) is rewritten as: 


q dA f +dA n = °1 S °3 d6 n dz n d0 f dz f ’ 


(A-9) 


where fti, ^2 and ^3 are given by: 


^1 = R n R f > 


( A-l 0 ) 


°* ■ 


[R f oos(e t - 8 n ) - R n ][R f - R n oos(0 f - 9 n )] 
[ R ? - 2 R n R f eos(0 f - 6 n ) + R 2 + (z f - z n ) 2 ] 2 


(A-ll ) 


= e f T f " e n T n * 


( A-l 2 ) 


Now, suppose that A^ is the finite furnace area "seen” by dA n . As shown in Figure A-2, 
Af consists of two separate regions on the outer cylinder plus the vertical edge of the 
heat barrier. Depending on the location of dA n , the area Af may also include one of the 
horizontal surfaces of the barrier. The outer cylinder portions of Af consist of a 
collection of elemental areas similar to dAf of Figure A-l. Thus, the radiative exchange 
between this portion of Af and dA n is described by equations (A-4) through (A-12). The 
heat barrier portions of Af, on the other hand, consist of elemental areas having variable 
Rf and different orientations with respect to dA n . Thus, Equations (A-4) and (A-7) should 
be rewritten for these elemental areas on the heat barrier. This can be avoided by 
realizing that the radiation exchange between dA n and dA^ (see Figure A-3) is equal to 


A-3 


that between dA n and dAf if dAf and dAf have the same thermal properties. That is 
because dA p "sees" dAf and dAf through the same solid angle. Thus, the horizontal 
surface of the barrier can be replaced by the outer cylinder region lying between points 
A and B in Figure A-3. It must be kept in mind that each element dAf in this region must 
be given the same thermal properties as its dAf counterpart. In a similar manner, the 
vertical edge of the heat barrier can be replaced by the region BC on the outer cylinder. 

The net energy exchange, per unit time, between a finite area, A n , on the inner cylinder 
and the furnace area, Af> seen by it can now be obtained by integrating Equation (A-9) 
over A n and A f : 


V-A = I I 0 1 n 2 a 3 d9 n dz n d9 f dz f ' (A - 13) 

f " A, A 
f n 

A positive result means that the energy leaves A f and reaches A n , while a negative value 
indicates a net energy loss from A n to Af. The area A n is described by the 
ranges (0 n 0 n 2 ) and (z n z R 2 ). Also, Af is defined by 
Of j-*- 0f 2 ) and (Zf z^ 2 ). Since the boundary nodes are arranged on the 
alloy’s surface with a 0 interval of 1 rad and a small vertical interval, Az n , it is 
convenient to write: 


9 n,l = 9 - <A-14> 

0 n,2 * 1 > < A - 15 > 

z n , 2 - z n,l = Az n ‘ (A-16) 

Also, as seen in Figure A-4, 0f j and 0f 2 are defined for each value of © n as 
follows: 


A-4 
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Fig. A-2 : Furnace View 


Fig. A-3 : 


Radiation Solid 
Angles 


(A-17 ) 


0 f,l = 6 n + sin _1 (R n /R f ) - w/2, 
0 f,2 = 6 n" sin _1 (R n /R f ) + n/2. 


( A-l 8 ) 


Finally, as suggested by Figure A-3, j= -°° and 2 = +0 °. Since ft^ is constant 
and ft g changes only with vertical distance, Equation (A-13) is rearranged to: 


z f , 2 z n, 2 ®n,2 0 f,2 


q A,-»-A_ _ fi l / 


l f ~n 


/ *3 J / 


ft 2 d0j d© n dz n dzp 


z f,l z n, 1 0 n,l 0 f,l 


(A-19) 


The result of the innermost integration, as suggested by Figure A-4, cannot be a function 

of 0_ . Also, since 0„ 0 - 0„ , =1, and Az is small such that: 
n n , 2 n , l ’ n 


z n , 2 

j $ dz n = Az n $ , 
Z n , 1 


( A-2 0 ) 


it follows that: 


z f , 2 0 f , 2 


q A = ft, Az n j ftj j ft„ d0, dz 

z „ , 0 


W 1 


3 J “'2 f 

f,l 6 f,l 


( A-2 1 ) 


Using the transformation ®f _0 n =T1 » the innermost integral is transformed to the 
following form: 


f,2 

J 

'f , i 


n 2 de : 


' / 


a cos n + b costi + c 


dTi = r(z n ,z f ), 


n j (d cosq + e) 


( A-22 ) 


A-6 



where: 


a = c = d/2 = -R n Rf , 
b = + R 2 , 

e = R n * R f * (z f‘ z n )2 
hj= sin *(R n /Rj) ” ” |r / 2 

n 2 = -sin _1 (R n /R f ) + tt/2 
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( A-2 3 ) 
( A-24 ) 
( A-2 5 ) 
( A-26 ) 
( A-27 ) 


It is clear that, for constant R n and Rf, the integral (A-22), denoted by r ( z , z ^ ) , 
is a function of the axial locations Zf and z n only. Taking R n =5 and Rf = 17.5, the 
integral (A-22) is solved numerically using the trapizoidal rule and the results are plotted 
in Figure A-5 for different values of z n and Zf. It is found that all the curves in Figure 
A-5 can be accurately duplicated by the single function aj sech (a 2 x) where x is equal to 
Zf - z n . The constants aj and a 2 are, of course, functions of R n and Rf. For the values 
of R n and Rf given above, aj and a 2 are found to be 0.0086 and 0.153, respectively. 
Equation (A-21) can now be rewritten as: 

z f , 2 

V-A * a l R 1 iz n / a 3 

f n Z 

Z f,l 

The above integral must be evaluated numerically for each A R which denotes a 
cylindrical surface area associated with a certain boundary node. By chosing R n , the 
area A n can be positioned anywhere between the alloy’s surface and the glass-air 
interface. In Equation (A-28), z n is the axial location of the centroid of A n which 
coinsides with the axial location of the associated node. The limits z^ j and Zf 2 can be 


sech a 2 (z£- z n ) dz^ . (A-28) 


A-7 





taken as large negative and large positive values respectively instead of -® and +<*> . 
This is possible because r ( z fi , z ^ ) diminishes at a fast rate as 
Zf - z n increases. 

The integral (A-28) can be evaluated numerically by dividing the region 
(z^ i -*• z^ into a larg® number of elements with height Az^ and centers at a 
variable Zf. The integral (A-28) takes the following discrete form: 


q A.H.A = a l 
f n 


ft, Az„ Az 
1 n 


Y sech[a 2 (z f - z n >] . 

z f 


( A-29 ) 


Finally, the radiation flux in the vicinity of the node is obtained by dividing Q by 

A f n 

the area A n : 


Q r ” q A,-»-A 1 A n * 
f n 


( A-30 ) 


By positioning A n at the middle of the glass wall, the above Q r becomes representative 
of an average heat flux exchanged between the furnace and the ampule by radiation. 
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